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I. INTRODUCTION 



II. GROUND STATE OF THE HUBBARD MODEL 



A short-range order, i.e. strong short-range correla- 
tions, is an intrinsic feature of strongly correlated Fermi 
systems. It was observed in a metal phase of V2O3 - 
a famous Mott- Hubbard system high-temperature 
superconductors (2), and heavy fermion systems H-f§. 
Principal aspects of the short-range order can be investi- 
gated within the Hubbard model The exact solu- 
tion of the Hubbard Hamiltonian is known for ID chain 
|^o| . Great simplifications appear in infinite dimensions 
because spatial correlations do not play an important 
role in this limit fjj]| . As for lattices of intermediate di- 
mensions, which are of great practical importance, one 
has to follow various numerical and analytical approx- 
imations [ p2||l4| ]. Recently a new variational approach 
to the ground state of the Hubbard model was pro- 
posed |T^]. In addition to intrasite correlations the trial 
wave function of Gutzwiller's type |gj contains nearest- 
neighbor correlations in an explicit form. In contrast to 
Rcf. 14|, Kikuchi's (cluster variation) method ]lq] was 
used to evaluate the ground state energy. It was shown 
that the short-range correlations affect significantly the 
ground state of the Hubbard model. A comparison of this 
result with the variational Monte Carlo method (VMC) 
jl7j based on the Gutzwiller trial wave function and 1 / D- 
expansion in the dynamical mean-field theory p8| shows 
that the latters underestimate significantly the ground 
state energy at intermediate coupling. 

A heavy fermion behavior usually arises from an inter- 
play between a lattice of localized /-electrons and itin- 
erant electrons. In the Kondo limit, such a system be- 
comes a Kondo lattice p9| . In many cases the itiner- 
ant subsystem of the Kondo lattice is formed by a nar- 
row d-band where the short-range Coulomb interaction 
between itinerant electrons is considerable. There is a 
lot of examples of this kind fL9|j . Therefore it is reason- 
able to describe the itinerant subsystem by means of the 
Hubbard Hamiltonian. Thus, we come to the Kondo- 
Hubbard lattice model [Q. It was shown by means of 
neutron scattering experiments that strong short-range 
antifcrromagnetic (AFM) correlations exist in Kondo and 
Kondo-Hubbard lattices ||-^|. They play an important 
role in the heavy fermions behavior. In this paper, I apply 
the variational theory of Ref. [[l6| to the Kondo-Hubbard 
lattice model. In Sec. II, the variational technique is in- 
troduced. The ground state energy of the paramagnetic 
(PM) phase of the Hubbard model at half-band filling is 
calculated. In Sec. Ill, the technique is generalized to the 
Kondo-Hubbard lattice model. The ground state energy, 
correlation functions, and effective mass are calculated 
at half-band filling for (i) spin 1/2 and (ii) spin 5/2 (/- 
spin) Kondo-Hubbard lattices. The results are discussed 
in Sec.IV. 



A. Trial Wave Function 

Let us consider a lattice with one orbital per site 
and restrict ourselves to nearest-neighbor hopping only. 
Then, the Hubbard model has the following form 



(1) 



(ij), 



where a irT {a,j a ) is the creation (annihilation) operator of 
a fermion of spin a =\,\ on the i-th lattice site, (ij) 
denotes a pair of adjacent sites, rn a = a} ia aj a . 

The Gutzwiller trial wave function Q gives a good ba- 
sis for a variational analysis of the Hubbard model ground 
state in infinite dimensions. This trial function can be 
written in the following symbolic form |21[| 



(2) 



where X — n^rt^, go is the real parameter taking a 
value in the interval [0, 1] if U > 0, |v?o) is the iV-particle 
wave function of non-interacting fermions, for instance 



m) 
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(3) 



keVj 



Here, aj^ denotes the creation operator of the Bloch 
state, k is the wave vector, and "vV CT is the space within 
the Fermi surface. 

To control nearest-neighbor correlations between 
fermions one can extend the Gutzwiller trial wave func- 
tion as |16l 



k'} = lH A i^> 



(4) 



where Pa are the projection operators onto all feasible 
configurations of a single site and a pair of nearest- 
neighbor sites, g\ are the nonnegative real parameters. 
In the PM phase, there are 4 such operators for intrasite 
configurations 

i 
i 
i 
i 

and 10 operators for nearest-neighbor configurations 
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All the operator Y\ and corresponding pair configura- 
tions are shown in Table |. 

From now on, we shall consider only lattices for which 
the total number of nearest-neighbors pairs is equal to 
zL/2, where z is the number of the nearest neighbors 
of a site and L is the total number of sites. Let us de- 
note normalized eigenvalues of the projection operators 
as x x |$) = L~ 1 X\ |$), y x |$) = (zL/2)^ 1 Y\ |$). The 
eigenvalues turn out to be related to each other by nor- 
malization conditions [122] 



= 1, ^2f3\V\ = 1 



(5) 



A A 

and self-consistency conditions |22| 

Vl +V3 + V4 + V5 =Xl, 

y-i + y-i + ys + yg^ x 4 , 
y4 + ye + y7 + y& = x 2 , 
V5 + y7 + y9 + yio = x 3 - 



(6) 



As concentrations of fermions of each spins are fixed 
there are the only independent parameter x\ and 7 inde- 
pendent parameters y\. In the case of half band-filling, 
additional constrains appear 



y\ = 2/2, ye = yw, 2/4 = 2/5 = ys = 2/9 



(7) 



which reduce the number of the independent parameters 
yx to 3. Assume that x\ = X4 = x, j/3, j/4, and yj are the 
independent parameters. Then, we obtain the final form 
of the generalized trial wave function of the PM phase at 
half band-filling 
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(8) 



To elucidate the physical meaning of the trial func- 
tion (@) let us rewrite the initial wave function of non- 
interacting fermions as a superposition of configurations 



|^ > = ^^r|r), 

r 

|r> = n«tl°) 

where Ar is the complex amplitude of the configuration 
r, |0) is the vacuum state. Then, we obtain 



M=Y,Vo ro 9 2 3 Dr3 9 8 4 D ™9 2 7 D "A r \T) 
r 



(9) 



Here -Dro denotes the number of doubly occupied sites in 
the configuration T, -Dr3(4.7) are the numbers of configu- 
rations of nearest-neighbor pairs corresponding to the op- 
erator I3, Y4, and Yj in the configuration T. Since other 



operators Yx are related to them by Eqs.(||), (||), and (Q), 
all the possible configurations of nearest-neighbor pairs 
are taken into account explicitly. 

Since the operator F = g^g^ 3 3 <? 4 g 7 7 is a poly- 
nomial of riia, it commutes with an operator of parti- 
cle alternation. Thus, the trial wave function (^) is an- 
tisymmetric. F is also invariant under the operations 
transforming lattice into itself, namely, translations, re- 
flections, rotations and inversion. This means that all 
the symmetries of the initial wave function remain valid 
for the trial wave function. 



B. Ground State Energy of the PM phase 

It is well known that correlations between fermions 
of the same spin exist even in the initial state (^), i.e. 
at U = 0. That is why, first we should evaluate its 
norm using an equiprobable state. The norm of a cor- 
related state is calculated also using the equiprobable 
state. Finally, all expressions should be normalized to 
the norm of the initial state @ . Since operators F make 
up an operator group {F(g ,g3,g4,g 7 )F(g' ,g' 3 ,g' 4 ,g' 7 ) = 

F (9o9o, 939*3, 9494,9797) and FF^ 1 = 1 when 9\ ar e finite 
and nonzero) this procedure is equivalent to applying the 
Fermi sea state (|J) as the initial one in Eq.jft. Thus, the 
norm of any state generated by Eq.(||) is ju| 

= J2 W {x , ya ^ y7} g^ g r Ly3 9l zLyi 9 2 7 zLy7 
{x,y3,y4,V7} 

= ^{^3/3,1/4,2/7} ■ (10) 

{®,y3,yi,V7} 

A factor which is inessential for the further calculations 
is dropped in Eq. ( |l0| ) . The summation is performed over 
all the sets {x, y 3 , j/4, y 7 }. A lot of configurations are re- 
lated to the same set of the independent variables. Then, 
Wfx, 2/3,1/4,3/7} is the number of the configurations corre- 
sponding to the fixed set {x, y 3 , j/4, yj} or the weight of 
this set. To calculate this quantity we use Kikuchi's 
pseudo-ensemble method Jlj|. It should be mentioned 
that this method is practically exact for Bethe lattices 
and approximate for lattices with closed paths [ p2[ . Ac- 
cording to Kikuchi's hypothesis the weight of a set can 
be expressed as a product [pl| : 



W = TQ 



(11) 



Here, lower indices are omitted for the sake of simplicity. 
The quantity Q determines the number of arrangements 
of ten indistinguishable elements corresponding to Yx at 
zL/2 pairs, i.e. the multinomial coefficients 



(12) 



and 



3 



L\Y[ x (x x zL)\ 
{zL)\Y[ x (x x L)\ 



(13) 



is the fraction of proper arrangements in the pseudo- 
ensemble |22] ]. In Eqs.(p^) and ( fi3| ) the dependent vari- 
ables should be expressed in terms of x, 2/3, y 4l 2/7 as fol- 
lows 



x 2 = x 3 = 1/2 — x, 

y\=V2=x-y3- 2j/ 4 , 

ye = 2/io = 1/2 - x - y 7 - 2y 4 . 



(14) 



From now on, we use 2/2 and ye as brief notations of the 
expressions (|l4[). 

In the thermodynamic limit (L — > 00) we can retain, 
in the usual fashion [|l5U22| 1, only the terms of the se- 
ries © which are very close to the largest one that is 
the following condition should be valid {x, 2/3, j/4, 2/7} — > 
{x, j/3, 2/4, yi^MAX- All the other terms are exponentially 
small. Since R is a nonnegative function, it is convenient 
to search the global maximum of its logarithm rather 
than of itself. Let us transform the factorials involved in 
R by means of the asymptotic Stirling formula. Then, let 
us find the logarithm of R and retain the leading term 
on L only. A straightforward calculation yields 



1 In W = 2(z - 1) [x In x + (1/2 - x) ln(l/2 - x) 



~z(y 2 In 2/2 + 2/3 In 2/3 
+2/6 hi 2/6 + 2/7 hi 2/7)- 



42/4 In 2/4 



(15) 



The domain of function ( jl5| ) is limited by conditions (|5|) 
and (|). It can be shown that the gradient of the func- 
tion at boundaries is directed inwards to the domain. 
Therefore the global maximum of R is an internal one 
and conditions 



djhiR) 
di]\ 



0, 



(16) 



where rj\ — x, 2/3, J/4, 2/7: ar e necessary for the global max- 
imum. They lead to the following system of equations 
that relate gi to x and yi 



90 = 


( 1/2 rr 


( x-y 3 -2y 4 \ 2/2 




\i/2 - s -2/7- 22/4 y 


9l = 


2/3 




x - 2/3 - 2?/4 ' 




9t = 




2/1 


(1/2 - z - 2/7 - 


- 2y 4 ) (x - 2/3 - 22/4) ' 


.97^ 


2/7 




l/2-x-y7- 


22/4' 



(17) 



It should be mentioned that L _1 lni? is rigorously a 
convex upwards function of 2/3, 2/4 , 2/7 at any fixed x. This 
means that, in effect, we search the maximum of a func- 
tion of the only inexplicit variable. 

To calculate the ground state energy of the Hamilto- 
nian (|l|) we need to evaluate the density matrix of the 
first order using the trial function (||) 



P 



1 

T 



E (<4v a 3° 

>4>>.rT ^ 



H.c 



<ij>,a 



WW) 



(18) 



Here, there is a significant complication as compared to 
the Gutzwiller trial wave function. While a fermion hops 
from site i to site j, the initial configurations of pair i — j 
and adjacent pairs (e.g. i — k, j — I) change (see Fig. la). 
Let us fix a configuration of pair i—j and adjacent pairs 
i — k,j — l and calculate function W of residual lattice by 
means of Eqs.(0), (H), and (H). The result is denoted 
by W. Then, a fraction of configurations containing the 
fixed fragment can be written as follows 



W -pr 

- W = ym{{ 



y{ki) 

x i0 



nfe) 



(19) 



where yu^ means some y\ corresponding to pair (ij). 
The term of the density matrix which comes from the 
transition from configuration 1 to configuration 2 takes 
the form 



riff*(2) w\i) 

Xli9i{l) W 



(20) 



where the first factor is the ratio between amplitudes of 
configuration 1 and 2, i.e. </i(l) and <?i(2) correspond to 
the configurations 1 and 2. In general, the procedure is 
similar to the Gutzwiller method || but the only param- 
eter go enters into Eq .(p0[ ) in the latter case. 

By means of Eqs. (|19|) and (^TJ) one can sum up over 
all the configurations and calculate the density matrix 



2y 4 (aia 2 ) 2 



2/357 2(*-l) 
a l 

9o93 



2/72/053 2(z-l) 

a 9 

97 



(21) 



where 



and 



«i 



a 2 



2/22/4 + 2/32/4.93 



2/4 (97 + 1) .94 1 



2/654 + 2/75457 



+ 2/4 (53 + 1)54 1 



1/2- 



The first term of Eq.(|2l|) describes transitions which 
do not change the total number of doubly occupied sites. 
The second and third terms are due to transitions corre- 
sponding to annihilation or creation of a doubly occupied 
site. Let us exclude parameters 5, by means of Eqs.(fi"7|). 
Then, after straightforward simplifications we obtain 



P = 8(2/4 + V 2/32/7) 

X (V2/2 + VV3 



x(l/2 



Y7)r 



(22) 



Finally it is convenient to present the total energy of 
Fermi system in Gutzwiller 's form 
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i (^W) 



(23) 



where </ = p° is the density matrix (22) at U = 



so = 2V F i / e k rfk 

is the average energy of the non-interacting fermions. 
First we calculate p° by minimization of the ground 
state energy at U = 0, i.e. p° = mm{ x ^ yim } (p). 
The ground state energy at nonzero U is determined 
as W3iit x y3>y4i Vr \ (E). The function E turns out to be 
smooth, without singular points within the domain of 
the function and its minimum is easily found numeri- 
cally. I used a refined Nelder-Mead simplex algorithm 
for the minimum search. 

The ground state energy of the PM phase calculated 
by this method is shown in Fig. 2: (a) a one-dimensional 
chain (z = 2) with the dispersion law = —2cosk xl (b) 
a square lattice (z — 4), — — 2[cos/ca; + cosk y ], (c) a 
simple cubic lattice (z = 6), £k = — 2[cosfc;2; + cosk y + 
coak z ]. Symmetric and antisymmetric correlation func- 
tions of the nearest neighbors 

G s = (n T n T )' + (njnj.)' = 2 (y 2 + 2y 4 + y 6 ) , 

G a = {n ]ni Y + ( ni n r )' = 2(y 2 + 2y 4 + y 7 ) (24) 

are shown in Fig. 3 for the same lattices as in Fig. 2. The 
prime in Eqs.(|24|) denotes the averaging over nearest- 
neighbor pairs only. Further details of the calculations 
can be found in Ref. U&. This technique was also ap- 



plied to the AFM phase of the Hubbard model 16 



C. Low-energy spectrum and effective mass 

To investigate the energy spectrum of low-lying exci- 
tations in a model with the total energy of the form (23), 
one can use the Fermi liquid approach of Ref. plj]. Let 



Eg be the ground state of the system (B3J). Then, let us 



create a new initial state \ifk 
trial wave function 



= a^ pa .a p i a \<po) and a new 



(25) 



where k = p p', p and p' are wave vectors lying above 
and below the Fermi surface, correspondingly. \tpo) is 
the ground state of non-interacting fermions. It should 
be noted that the new trial wave function has the same 
number of fermions of each spin as the initial wave func- 
tion (the canonical ensemble) . Since the operator on the 
right hand side of Eq. (||) is translationally invariant the 
new trial wave function has the well-defined wave vector 
k. Then, we perform the procedure developed above to 
determine the minimum energy E g + 8E-^ a corresponding 
to the excited state |Vw), i.e. we find a new stationary 
solution of the Hamiltonian ([!]). It is easy to see that 



SE^cr is small, of the order of 1/N. We express the en- 
ergy variation as 



SEto = q 6E^ 



-£o 



dq 
dx 



Sx 



E 



dXi 



U Sx 



(26) 



where = 2/3,2/4, y-; and 5E^ a is the energy variation of 
non-interacting fermions corresponding to the excitation 
\ifka-). Since the ground state is minimal the following 
conditions are valid 



£07^ + ^ = 0, 
ox 



dq 



0. 



(27) 



Combining Eqs.(|2q) and ( |27| ) we find that SE^ 
q 5E^ a for any trial wave function generated by Eq.(^5| 
i.e. the low-lying energy spectrum of our model is 



£k<T 



9 4. 



where e ko . is the energy spectrum of non-interacting 
fermions. Retaining the terms of the second order in- 
finitesimal in the expansion ( |26j ) we obtain the effec- 
tive Fermi liquid theory One can perform the 
similar calculation with the grand canonical ensemble 
(Iv'kcr) = a ko-l < ^o)) but it would be more complicated be- 
cause, in this dependence of q on the number 
of fermions has to be taken into account. The effective 
mass at the Fermi surface is m = q^ 1 where the effec- 
tive mass of non-interacting fermions is assumed to be 
too = 1. This result is in an agreement with the phe- 
nological Brinkman-Rice approach p3[ and slave-boson 
treatment p4] . The effective mass in the half- filled Hub- 
bard model is shown in Fig. 4 as a function of U. A de- 
tailed discussion of the excitation spectrum will be pre- 
sented elsewhere. 



III. GROUND STATE OF THE 
KONDO-HUBBARD LATTICE 

The Kondo-Hubbard lattice model can be expressed in 
the following form [M 



Hkh — Hh + J ^ S -Sf 



2 \Si+S£ 



51 QC 



(28) 



where is the Hubbard Hamiltonian, Sf is the spin 
operator of an itinerant fermion at site i, S\ denotes the 
spin or the total angular moment operator (/-spin) de- 
pending on the nature of the localized state. We consider 
the PM phase at half band-filling. It is convenient to rep- 
resent the Kondo-Hubbard lattice as an equivalent lattice 



5 



in Fig. lb. Here, we obtain a new sort of nearest-neighbor 
pairs, namely, itinerant fermion - localized fermion on the 
same site (open circle - black circle in Fig. lb). It should 
be noted that there are no additional closed paths in 
the lattice as compared to Fig. la. That is why, we don't 
bring about additional simplifications as compared to the 
Hubbard model. The general form of the trial wave func- 
tion of the Kondo-Hubbard lattice can be presented as 



IV' 



KH 



9? r \4>h) \<pi) 



: n^ivc)iw) (29) 



where k'cfn) is the initial wave function of itinerant (c) 

and localized fermions (I), P\ are the projection opera- 
tors for the Hubbard model (Q), \tpu) is the trial wave 
function of the Hubbard model and 9? = — S l zi S zi is 
the new projection operators for itinerant fermion - lo- 
calized fermion pairs. Index z denotes the projection on 
z-axis. \tpi) is the PM phase without correlations (all spin 
configurations are equiprobable). In the next two subsec- 
tions, we shall define concretely the trial wave function 
( p9| ) for two cases. 

A. 5 = 1/2 Kondo-Hubbard lattice 

There are three eigenstates of operator S'S^ in case of 
the spin 1/2 localized state (singlet, triplet and Sf z = 0, 
see Table Let us introduce the eigenvalues of the 
operator ro, r\ and r 2 - The self-consistency and nor- 
malization conditions (||) , (|^) remain valid and new ones 
appear 



r = x, n+r 2 = 1/2 - x. 



(30) 



From this it follows that a new independent parame- 
ter r = t*i appears in addition to the set describing the 
Hubbard subsystem (x, 7/3, 7/4, 2/7). Then, the trial wave 
function takes the form 



\1pKn) = gf \iPh) \<pi) 

5R X 03% B4Y4 3 7 Y 7 1 \ 1 \ 
= 9r9o 93 94 97 \fc) \<Pl) 



(31) 



and its norm is 



<v#) = E 

{r,x,y 3 ,y 4 ,y 7 } 



W, K X W, H x 

{r,x\ {x,y 3 ,y4,,y7i 



v r ALr 2Lx 2zLy 3 izhy A 2zLy 7 
x 9r 9o 93 9 A 9i 



E 



R 



{r,x,ys,y4,y 7 } 



(32) 



{r,x, 3/3,3/4,2/7} 



where W K ( S ' is the Kondo (Hubbard) weight of the set. 
W H is from Eqs.(£§) and @. The Kondo weight can 
be easily calculated 



W 



K 



[L(l/2-x)]\ 



(2xL)\ 

[(^f \[rL]\[L(l/2-x-r)}\ 



(33) 



Using Eqs. ( |32| ) and (33) we search the global max- 
imum of the norm following the procedure described in 
the previous section. We obtain the system of equations 
( |l6| ) where rj\ — r, x, 7/3, 7/4, j/7. The equations for param- 
eters (73, <?4 and gj remain the same as in system ^Ll\). 
For the other parameters we obtain 



5o 



1/2 - x 



2y/r{l/2-x-r) 



1/2 - x 



x-Vz- 2y 4 
1/2 - x - yr - 2y 4 



z/2 



1/2- 



(34) 



The total energy of the Kondo-Hubbard lattice in- 
cludes Hubbard and exchange parts 

E = - -TTTT\ = q£o + xU + J (p zz + p± ) (35) 

where q = p#/ p H , p° H is the Hubbard density matrix at 
U = 0, p zz , p± are the density matrices corresponding to 
zz and spin-flip interactions. Since there are new bonds 
in the lattice in Fig. lb, the Hubbard density matrix is 
different from that of the Hubbard model. Thus, instead 
of Eq.@ we get 



PH 



where 



2y 4 (a 1 a 2 ) z hb 2 



V79093 2(z-l),2 



.97 



(36) 



h =rg r 1 +g r (l/2-x-r), 

b 2 = 2 (ffr +9r X ) ■ 

Straightforward calculations of exchange terms give 
Pzz = ~(2r + a;-l/2), 



(37) 



P± 



2j/4 + V697 + V79 7 

(1/2 -xf 



Expressions (|35|), ( p6| ) and ( p7| ) present the total en- 
ergy of the Kondo-Hubbard lattice in an analytic form 
as a function of independent variational parameters r, x, 
V3i Hi > V7- The ground state energy is the global mini- 
mum of Eq.(35) with respect to these parameters. The 
spin nearest-neighbor correlation function of the itiner- 
ant subsystem 



G B = 



(S c zi S c zj )' = (G s - G a ) /4 = ~ (t/6 - y 7 ) 



(38) 



and the spin correlation function of localized nearest 
neighbors are shown in Fig. 5. The last is calculated by 
means of the superposition hypothesis (22| 



G 



2x - 4r) (y 6 - y 7 ) . 



(39) 



Following the expansion (g6|) we determine the effective 
mass of the itinerant fermions asm = q 1 - It is plotted 
against J in Fig. 6. We have also plotted the fraction of 
lattice sites occurred to be in the Kondo singlet state (2r) 
in Fig. 7. 



B. S = 5/2 Kondo-Hubbard lattice 

Bearing in mind cerium Kondo-Hubbard lattices 
(Ce 3+ ) we generalize the technique developed above to 
spin 5/2 (the total angular moment or /-spin]. The crys- 
tal field is neglected in the Hamiltonian (EST). Crystal 
field effects will be discussed briefly in the next section. 
Let us compile a table of all intrasite configurations (see 
Table H|). Since the Hamiltonian ( p8| ) is rotation invari- 
ant, possibilities of configurations with Sf z — are equal 
each other, r$ = x/3. There are 6 configurations with 
Sf z =/= 0. They are bound by the normalization condition 



1/2 



(40) 



It follows that one of u is dependent (we take r§ as the 
dependent parameter). The general trial wave function 
( pl| ) remain valid for the spin 5/2. Then, its norm is 



E 

{ri,x,y 3 ,y4,y 7 } 



w< K X W, H x 



v n 4L Eli"-. (7-2i) 2Lx2zLy 3 SzL yi 2zLy 7 

x 9r 9o 93 9i 9i 



(41) 



In case of S = 5/2, the Kondo part of the configuration 
weight is reduced to 



W 



K 



(2xL)\ ([L(l/2-x)]\ 



l(xL)\} 6 I IL (nL)l 



(42) 



The necessary condition of the global maximum of R 
(plf), where tjx = r,,a;, 2/3,1/4,2/7, leads to a system of 
nine equations. In particular, there are five equations for 
independent 



2(12-2i) 
It 



(43) 



where r§ = 1/2— x— rj. From here one can see that, in 
effect, there is the only independent parameter r. Let it 
be r = n . This is no wonder because the only variational 
parameter r enters into the trial wave function. The five 
equations for ri gives the nonlinear equation relating r 
and g r 



1 



9r 



9r 



- 9r 



-12 



9r 



-16 



9r 2 °) 



1/2- 



(44) 



The equation for go transforms to 



9o 



1/2- 



0>Jrr§ 



1/2- 



x 



x 

2/3 - 2j/ 4 



1/2- 



x ■ 



2/7 - 22/4 



r/2 



(45) 



The equations for 2/3, 2/4, 2/7 are t ne same as for the Hub- 
bard model @. 

The density matrix falls into the Hubbard, zz and spin- 
flip terms in the same manner as for the localized spin 
1/2. The general form of the Hubbard term is Eq.(|3^) 
where 

2—1 



1/2 



Straightforward calculations give the zz and spin-flip 
terms 



P± 



1 6 

i=i 
-(1 - 2x 



r 9 ; 2 °)9r 2 - 



(46) 



The total energy is written in the form of Eq. ( |35| ) . The 
ground state energy is determined by the minimization 
over the independent variational parameters. There is a 
complication of the numerical procedure because in the 
present case the variational parameter g r can't be ex- 
pressed analytically in terms of r. In other respects the 
calculations are similar to that discussed in previous sec- 
tions. The spin correlation function of itinerant fermions 
( |38| ) and that of localized states 

Gi = ( x S zi S Z j') [(S^jS^j)] 
- 1 (2/6 - 2/4) r 2 



4 

x(5 [1 



-201 



■3[ 5 ,7 4 



9, 



-161 



9 r 



9r 



12\ 



(47) 



are shown in the insert in Fig. 5. The effective mass and 
the function 2r are plotted in the inserts in Fig. 6 and 
Fig.7. 



IV. DISCUSSION AND CONCLUSIONS 

The trial wave functions used in the present approach 
(^), (^9|) have a remarkable property. The operator on 
the right side of Eq.([|) commutes with the operators of 
the crystal point group (rotations, reflections, inversion) 
and the translation group. That is why the trial wave 
function retains all the symmetries of the initial wave 
function. In contrast to the Gutzwiller wave function, we 
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have incorporated the nearest-neighbor correlations into 
the trail wave function. That is why, a local structure sur- 
rounding an atom appears. A similar trial wave function 
was used in Ref. The ground state energy was 

evaluated there by means of an expansion of exponent 
exp(— "fY\). To perform this a systematic diagram repre- 
sentation was introduced This approach turned 
out to be very successful in quantum chemistry for small 
molecules because one obtains an excellent results taking 
into account the first terms of the series. At the same 
time, the eigenvalues of Y\ are approximately propor- 
tional to the number of lattice sites. That is why, it is 
hardly possible to use this approach in the thermody- 
namic limit (L — > oo) for the strong short-range order. 
Let us consider a set of approximations in the framework 
of Kikuchi's method, namely the Gutzwiller approxima- 
tion (a cluster consists of the only site) as the first one, 
a pair Kikuchi's approximation (a cluster consists of two 
sites) as the second one, and etc. On each step, a cluster 
becomes larger. It was rigorously proved for 2D square 
lattice and 3D cubic lattice that in the thermodynamic 
limit we approach the rigorous solution while a cluster 
size increases [^7j. This allows us to go beyond the infi- 
nite dimension limit. 

We calculated the ground state energy of the Hub- 
bard model in the PM phase for a one-dimensional chain, 
square and simple cubic lattices (see Fig. 2). The result 
for a one-dimensional chain is very close to that obtained 
by an analytic investigation of Gutzwiller wave function 
p6fl . Let us mention that the method presented above 
gives much lower the ground state energy in the AFM 
phase The results for the square and simple cu- 

bic lattices are compared with that of the VMC method 
in Fig.2b,c. Since the VMC method is based on the 
Gutzwiller trial wave function the difference between this 
approach and our method is due to effects of the short- 
range correlations. It can be seen that near Uc = 8 |eo| 
(i.e. the critical value of U in the Gutzwiller approxima- 
tion) the ground state energy of the trial wave function 
Eq. (H) is substantially lower (two or tree times) than that 
obtained by the VMC method, i.e. the short-range order 
considerably reduces the ground state energy. This re- 
sult is in contrast to the 1 /D-expansion in the dynamical 
mean field theory |l8j and shows that the perturbation 
theory methods are hardly suitable here. 

The narrow quasi-particle band (see Fig. 4) is appeared 
at the Fermi level similarly to the results of Gutzwillcr's 
approach Pj2l]|. At large U, the effective mass become 
linear on U. The slope of the function m(U) increases 
while the lattice dimension becomes larger. The limit of 
infinite dimensions of the ground state ([)]) was investi- 
gated in Ref. |l6|. While the lattice dimension increases, 
the ground state energy approaches to Gutzwiller's so- 
lution. We also note that an exchange hole exists in 
the variational solution even at U = (G s < 0.5). 
It increases monotonically while U grows (exchange- 
correlation hole). Unlike the well-known Hubbard III 



solution the short-range AFM correlations do not disap- 
pear in the limit \t\ /U <C I but tend to a certain con- 
stant value (see Fig. 3). In this limit, the Hubbard model 
at half band- filling reduces to the spin- 1/2 Hiesenberg 
model. The residual AFM correlations in the PM phase 
are consistent with the results of the ground state study- 
ing of the Hiesenberg model p8| . It was also shown in 
Ref. [[l6| that in the AFM phase, the both methods (the 
VMC method and the present one) give almost the same 
ground state energy, i.e. in the presence of the long-range 
order, the short-range correlations are inessential. 

In the Kondo-Hubbard model we observe two differ- 
ent types of behavior depending on value of J. At small 
J, the exchange term weakly affects the ground state en- 
ergy and effective mass. It should be noted that the AFM 
correlations between localized states at a pair of nearest- 
neighbors sites (Gi) increase with increasing of J. This 
effect appears because the short-range AFM correlations 
between band fermions exist at J = and at small J 
the correlations between localized states follow the cor- 
relations of band fermions. The growth of the exchange 
leads to suppression of the band correlations. In Fig. 5, it 
is seen a transition to another regime (large J). For sim- 
ple cubic lattice we have obtained a discontinuity of G c 
and a sharp turn of Gi for both S l = 1/2 and S = 5/2 
cases, i.e. a phase transition. A smooth crossover is ob- 
served in ID chain and square lattice at this point (see 
Fig. 5). At present, it is hardly possible to establish a 
kind of the phase transition. We can rule out the first 
order transition only because the ground state energy is 
smooth (the derivative of the ground state energy is con- 
tinuous). At higher J, the correlations between localized 
states start to decrease with increasing J. It should be 
mentioned that the behavior of the spin correlation func- 
tion Gi is similar to that ensued from Doniach's phase 
diagram despite the fact that the RKKY interaction is 
not included into the variational theory. 

The Coulomb interaction between itinerant fermions 
U influences the ground state in two ways. On the first 
hand, it reinforces the short-range order and increases 
the AFM correlations between nearest-neighbor localized 
states. On the other hand, at large U the effective quasi- 
particle band gets smaller that favors the Kondo regime. 

At large J, the ground state energy approaches asymp- 
totically the energy of the pair singlet state (|J for 
S l = 1/2). The probability of the Kondo singlet (2r) 
at a site goes to 1 (see Fig. 7) for S l = 1/2 and to some 
fixed value (independent on the lattice dimension) for 
S l = 5/2. The difference is due to the fact that, in the 
first case, the ground state is a superposition of two an- 
tiparallel states (-^(Ifl) — an d, in the case of a 

pair S l = 5/2, S° — 1/2, it involves all the possible 
states (\S l zl = f) \S c zi = -|>, \S l zl = §) \S c zi = -f) and 
etc.) for the total energy to be minimal. At large J 
the effective mass of itinerant fermions increases drasti- 
cally and tends to almost linear behavior with variation J 
(to oc a J where a is the constant depending on the type 
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of lattice) as can be seen in Fig. 6. This is the Kondo 
lattice regime. In the case of S l = 5/2, the effective mass 
is significantly larger then for S = 1/2. In a real system, 
m can't increase infinitely and a localization of fcrmions 
should occur due to disorder or temperature effects but 
this is beyond the scope of the present paper. 

The crystal field can be taken into account in the vari- 
ational scheme. It can be easily done in the limit of 
the strong crystal field when we consider the lowest level 
of the multiplet only. In the general case, we have to 
use additional variational parameters to control differ- 
ent populations at multiplet levels. For instance, if the 
ground state of Ce 3+ ion is splitted by the crystal field 
into a doublet and a quadruplet (e.g. T*? and T 8 states), 
we should apply the only additional variational parame- 
ter to the trail wave function. If it is splitted into three 
doublets, we have to use two additional parameters and 
etc. 

In the framework of the variational theory presented 
above the strongly correlated coherent metal state ap- 
pears in a natural way similarly to the Gutzwiller the- 
ory. This provides a good basis for investigations of 
strongly correlated metal systems and dense Kondo sys- 
tems by means of this theory. The main shortcoming 
of our approach is neglect of closed loops while we treat 
the correlations on a lattice. Nevertheless, it is consid- 
ered that Kikuchi's method yields a good approximation 
if the correlation length is not greater than a cluster size 
p9[ . The well-developed cluster variation method allows 
to include short closed loops, which are most important, 
into considerations |}0| . The short-range correlations are 
observed directly by the neutron scattering. The AFM 
correlations in the strongly correlated metals like V2O3 
and Kondo lattices turn out to be strong but very short 
|j])|r[§- That is why the closed loops are inessential in 
this substances. Let us mention that one can calculated 
the neutron cross-section and dynamical susceptibility 
from the results obtained above. At the same time, it 
is hardly possible to apply the present theory to sys- 
tems with lengthy AFM correlations (e.g. strongly un- 
derdoped high temperature superconductors Q). 
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VI. FIGURE CAPTIONS 

Fig.l. (a) A fragment of z — 4 lattice, (b) represen- 
tation of z — 4 Kondo-lattice. Itinerant and localized 
states are shown as light and black circles, correspond- 
ingly; double lines denote the exchange interaction. 

Fig. 2. The ground state energy of the Hubbard model 
at half-band filling, (a) A one-dimensional chain: the 
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Gutzwiller solution (1), present study (2), the exact so- 
lution; (b) square and (c) simple cubic lattices: the 
Gutzwiller solution (1), the VMC method (2), the present 
study (3). 

Fig. 3. The symmetric G s (dashed lines) and antisym- 
metric G a (solid lines) correlation functions of the Hub- 
bard model for a one-dimensional chain (1), square (2) 
and simple cubic lattices (3). 

Fig. 4. The effective mass in the half-filled Hubbard 
model: a one-dimensional chain (1), square (2) and sim- 
ple cubic lattices (3). The dashed lines are guides to eye. 

Fig. 5. The spin correlation functions G c (dashed lines) 
and Gi (solid lines) for the half-filed S — 1/2 Kondo- 
Hubbard model: (a)_a one-dimensional chain, (b) square 
and (c) simple cubic lattices; U — 0.5U C . The correla- 
tion functions for the half-filed S = 5/2 Kondo-Hubbard 
model are shown in the insert. 

Fig. 6. The effective mass in the S = 1/2 Kondo- 
Hubbard model: a one-dimensional chain( 1), square (2) 
and simple cubic lattices (3); U = 0.5U C . The dashed 
lines are guides to eye. The effective mass for the half- 
filed S — 5/2 Kondo-Hubbard model is shown in the 
insert. 

Fig. 7. The probability of the Kondo singlet state in the 
S =1/2 Kondo-Hubbard model (2r). The probability of 
\S l zi = ±5/2) \S c zi = Tl/2) states in the S = 5/2 Kondo- 
Hubbard model is shown in the insert. 

TABLE I. Pair projection operators, corresponding config- 
urations and the degeneracy factor 
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TABLE III. Intrasite configurations of the spin 5/2 
Kondo-Hubbard lattice 

Configuration Eigenvalue Degeneracy 



Itinerant state S z Localized state J z f3i 

any state x 2 

1/2 -5/2 n 2 

1/2 -3/2 r 2 2 

1/2 -1/2 r 3 2 

1/2 1/2 r 4 2 

1/2 3/2 r 5 2 

1/2 5/2 r e = 1/2 -x 2 



TABLE II. Intrasite configurations of the spin 1/2 
Kondo-Hubbard lattice 



Configuration 


Eigenvalue 


Degeneracy 


Itinerant state Localized state S z 
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